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METHOD liJUD APPARATUS FOR MULTIVARIATE SPACE PROCESSING 

BACKGROUND OF THE INVENTION 

5 1. Field of the Invention 

The present invention relates to a multivariate space 
processing technology and multivariate space processing method 
and apparatus, and more particularly relates to a method and 
apparatus for predicting, interpolating and displaying 
10 multivariate data. 

2. Description of the Related Art 

Scientific visualization is often used as a tool to 
allow a user to intuitively grasp or understand complex 

15 phenomena and large-scale data. For example, in order to 
visualize time-varying complex three-dimensional phenomena 
such as ocean currents or tornados, density and other 
attributes of each point in the three-dimensional space are 
assigned for each voxel and time-varied values are computed 

20 sequentially and displayed, so that a general structure and 
motion of the large-scale data can be approximated. Besides 
grasping a phenomenon as a whole, visualization is also used 
as a tool to gauge the effectiveness of a hypothetical 
modeling simulation in a simplified and convenient manner. 
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In most cases, scientific visualization is implemented 
to allow a phenomenon or model to be understood in a general 
and simplified way. However, in actuality, the creation of a 
5 scientific visualization will generally involve a series of 
relatively complex processes such as setting up complex 
equations and functions, tracing time-varying values of each 
point in the n-dimensional space, and converting this large 
amount of data to a display space. Thus, the computational 
10 load for even a basic scientific visualization can be quite 
large, so that, in many cases, the amount of time and effort 
required to complete the scientific visualization adversely 
affects the pace at which the overall research processes can 
be conducted. 

15 



SUMMARY OF THE INVENTION 

The present invention has been made in view of the 
20 foregoing circumstances and an object thereof is to provide a 
technique by which the time-varying behaviour of three- 
dimensional objects can be visualized or displayed with a 
reduced computational load. Another object of the present 
invention is to provide a technique by which the structure of 
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or changes in a multivariate space having a dimension greater 
than three can be visualized or displayed in a simple and 
efficient manner. 

Hereinafter, the term "variate" is used to describe each 
of a plurality of variable quantities in a multivariate space. 
A variate may also sometimes be called a parameter and these 
terms may generally be treated in the same manner as the term 
"dimension" as the case may be. Further, when a variate is 
referred to as a variate of an object, the variate may also be 
defined as an attribute of the object. 

An embodiment according to the present invention relates 
to a multivariate space processing method. The multivariate 
space processing method includes: degenerating multivariate 
data into three predetermined variates (referred to as, for 
example, x, y, t hereinafter) ; determining a reference variate 
(referred to as, for example, t hereinafter) to serve as a 
reference among the three variates; acquiring a first two- 
dimensional space formed by the remaining two variates when 
the reference variate takes a first value; acquiring a second 
two-dimensional space formed by the remaining two variates 
when the reference variate takes a second value; and ' computing 
a matching between the first two-dimensional space and the 
second two-dimensional space. 

In particular, the ^'degenerating" may be such that 

MN-70022 



certain variate(s) or dimension (s) are simply cut off or, in 
other cases, may be such that a specific value is assumed for 
certain variate(s) or dimensions. For example, in a case 
where a three-dimensional space is degenerated to a two- 
dimensional space, the variate of height may be set to some 
constant (z=const.). Here, it suffices that the degeneration 
is such that an n dimension space is degenerated to a 
dimension which is less than n. For example, a three- 
dimensional space may simply be converted to an arbitrary 
plane ax+by+cz+d=0 . Moreover, a three-dimensional space may 
be converted to a curved surface such as ax^+by+c=0 . In any 
event, reducing the number of variates is called degeneration. 

As an example, consider the case where the reference 
variate t is time. Now, the remaining two variates x and y 
can be used to define a first two-dimensional space, sometimes 
called the first image, at time t=tO and a second two- 
dimensional space, sometimes called the second image, at time 
t=tl. If an image matching is computed between the first and 
second images, the approximate behavior of the two variates x 
and y in the region t=[tO, tl] can be determined. In this 
case, it suffices to compute a matching between the two- 
dimensional images, and the computational load can be reduced 
by selecting a suitable algorithm. The "base technology" 
described later provides one such suitable algorithm and 
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contributes to producing a suitable visualization. 

With this embodiment, a complicated multivariate space 
or object can be converted into a simple model by selecting 
three variates . As an example, consider that there are a 
5 large number of factors (variates) that go into determining a 
stock price. Initially, three variates of the stock price are 
tentatively determined based on a rule of thumb or an 
arbitrary selection after a search algorithm has enumerated a 
number of possible variates. Next, a first image and a second 

10 image are generated and an intermediate image thereof at, for 
example, t=t2 is generated by matching and interpolation. 
Thereafter, this intermediate image (also referred to as a 
virtual intermediate image hereinafter) and the actual stock 
price at the time t=t2 (also referred to as an authentic 

15 intermediate image) are compared to determine the suitability 
of the selection of the three variates. Thus, if the virtual 
intermediate image is close to the authentic intermediate 
image, it is concluded that the three variates have high 
importance. This can be used or further tested by determining 

20 if it is possible to analyze or predict the stock price at 
. another time by using these three variates. Further, the 
virtual intermediate image and the authentic intermediate 
image may also be compared while the three variates are being 
changed, while the selection of the variate serving as the 



reference (also referred to as the reference variate) is being 
changed, and/or while the value of the reference variate is 
being changed. 

By adopting the above-described methods while reducing 
the number of parameters to three, a stock price model can be 
simultaneously simplified and optimized. The methods may be 
stated from a more general perspective as: viewing a 
multidimensional phenomenon or object from various angles 
while reducing parameters one at a time, conducting an 
inspection to find an angle at which a most characteristic 
image is obtained, and repeating a degeneration to view the 
object at that particular angle. Using a simple example, 
consider a case where a three-dimensional object is projected 
as a "shadow picture" on a plane. There will generally be an 
angle at which the original object can be relatively easily 
grasped or understood from the shadow picture. A projection 
operation using this angle corresponds to the degeneration. 

The above-described matching may be performed by 
regarding the first two-dimensional space and second two- 
dimensional space as a first image and a second image, 
respectively, and computing pixel by pixel based on 
correspondence between a critical point detected through a 
two-dimensional search on the first image and a critical point 
detected through a two-dimensional search on the second image. 
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The matching may further include: multiresolutionalizing 
the first image and the second image by respectively 
extracting the critical points; performing a pixel-by-pixel 
matching computation on the first image and the second image 
at the same multiresolution ' level; and acquiring a pixel-by- 
pixel correspondence relation at a finer level of resolution 
while inheriting a result of the pixel-by-pixel matching 
computation from a matching computation at a different 
multiresolution level. 

Another preferred embodiment according to the present 
invention relates also to a multivariate space processing 
method that includes: acquiring a first image and a second 
image by projecting three-dimensional data on a predetermined 
plane; and computing a matching between the first image and 
the second image. By implementing this method, after the 
three-dimensional data of, for example, a tornado are 
projected on a two-dimensional plane to create the first and 
second images, intermediate states in a region of time t=[tO, 
tl] can be visualized by interpolating the two-dimensional 
images for time t=tO and time t=tl. Computationally, once the 
matching is completed, there will not be much time needed for 
the interpolation. Moreover, since the interpolation can be 
performed at any point in the region by arbitrarily varying 
the time t, the computational load becomes extremely low 

MN-70022 



compared to conventional visualization, in which a three- 
dimensional computation is performed sequentially. Thus, this 
embodiment is effective particularly when one wants to quickly 
visualize a phenomenon or simulation in as simplified and 
convenient a manner as possible. 

Still another embodiment according to the present 
invention relates to a multivariate space processing apparatus 
that includes: a preprocessing unit which degenerates 
multivariate data into three predetermined variates; a 
conversion unit which determines a reference variate from 
among the three variates to serve as a reference, acquires, as 
a first image, a two-dimensional space formed by the remaining 
two variates when the reference variate takes a first value, 
and acquires, as a second image, a two-dimensional space 
formed by the remaining two variates when the reference 
variate takes a second value; and a matching processor which 
computes a matching between the first image and the second 
image. In this embodiment, the apparatus may further include 
an intermediate image generator which generates an 
intermediate image of the first image and the second image by 
performing an interpolation computation based on a result of 
the matching computation. 

Still another embodiment according to the present 
invention relates also to a multivariate space processing 
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apparatus that includes: a conversion unit which acquires a 
first image and a second image by projecting three-dimensional 
data on a predetermined plane; and a matching processor which 
computes a matching between the first image and the second 
image . 

In the embodiments herein, the matching and, in 
particular, the matching using critical points may be an 
application of the techniques (referred to as the "base 
technology" hereinafter) proposed in Japanese Patent 
No. 2927350 owned by the same assignee of the present patent 
application. 

It is to be noted that the base technology is not a 
prerequisite in the present invention. Moreover, it is also 
possible to have replacement or substitution of the above- 
described components, elements, functions or processes in part 
or whole as between method and apparatus or to add components, 
elements, functions or processes to method or apparatus, and 
it will also be understood that the components, elements, 
functions or processes may be implemented by a computer 
program and saved on a recording medium or the like and are 
all effective as and encompassed by the present invention. 

Moreover, this summary of the invention does not 
necessarily describe all necessary features so that the 
invention may also be sub-combination of these described 
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features . 

BRIEF DESCRIPTION OF THE DRAWINGS 

Fig. 1(a) is an image obtained as a result of the 
application of an averaging filter to a human facial image. 

Fig. 1 (b) is an image obtained as a result of the 
application of an averaging filter to another human facial 
image . 

Fig. 1(c) is an image of a human face at p'^'°' obtained 
in a preferred embodiment in the base technology. 

Fig. 1(d) is another image of a human face at p'^'°' 
obtained in a preferred embodiment in the base technology. 

Fig. 1(e) is an image of a human face at p<^'^^ obtained 
in a preferred embodiment in the base technology. 

Fig, 1(f) is another image of a human face at p'^'^' 
obtained in a preferred embodiment in the base technology. 

Fig. 1(g) is an image of a human face at p*^'^' obtained 
in a preferred embodiment in the base technology. 

Fig, 1(h) is another image of a human face at p^^'^^ 
obtained in a preferred embodiment in the base technology. 

Fig. l(i) is an image of a human face at p*^'^^ obtained 
in a preferred embodiment in the base technology. 
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Fig. l(j) is another image of a human face at 
obtained in a preferred embodiment in the base technology. 
Fig. 2(R) shows an original quadrilateral. 

Fig. 2(A) shows an inherited quadrilateral. 
Fig. 2(B) shows an inherited quadrilateral. 
Fig. 2(C) shows an inherited quadrilateral. 
Fig. 2(D) shows an inherited quadrilateral. 
Fig. 2(E) shows an inherited quadrilateral. 
Fig. 3 is a diagram showing the relationship between a 
source image and a destination image and that between the m-th 
level and the (m-l)th level, using a quadrilateral. 

Fig. 4 shows the relationship between a parameterri 
(represented by x-axis) and energy Cf (represented by y-axis) . 

Fig. 5(a) is a diagram illustrating determination of 
whether or not the mapping for a certain point satisfies the 
bijectivity condition through the outer product computation. 

Fig. 5(b) is a diagram illustrating determination of 
whether or not the mapping for a certain point satisfies the 
bijectivity condition through the outer product computation. 

Fig. 6 is a flowchart of the entire procedure of a 
preferred embodiment in the base technology. 

Fig. 7 is a flowchart showing the details of the process 
at SI in Fig. 6. 

Fig. 8 is a flowchart showing the details of the process 
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at SIO in Fig. 7. 

Fig. 9 is a diagram showing correspondence between 
partial images of the m-th and (m-l)th levels of resolution. 

Fig. 10 is a diagram showing source hierarchical images 
generated in the embodiment in the base technology. 

Fig. 11 is a flowchart of a preparation procedure for S2 
in Fig. 6. 

Fig. 12 is a flowchart showing the details of the 
process at S2 in Fig. 6. 

Fig. 13 is a diagram showing the way a submapping is 
determined at the 0-th level. 

Fig. 14 is a diagram showing the way a submapping is 
determined at the first level. 

Fig. 15 is a flowchart showing the details of the 
process at S21 in Fig. 12. 

Fig. 16 is a graph showing the behavior of energy 
corresponding to f^'^'S' (X = iAZ) which has been obtained for a 
certain f^"''^) while varying X. 

Fig. 17 is a diagram showing the behavior of energy Cj-"^ 
corresponding to f<"' (77 = /A7) (z = 0,1,...) which has been obtained 
while varying t|. 

Fig. 18 shows a structure of a multivariate space 
processing apparatus according to an embodiment of the 
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invention. 

Fig. 19 is a flowchart showing a multivariate space 
processing method according to an embodiment of the invention. 

Fig. 20 shows more detail of S300 shown in Fig. 19. 
5 Fig, 21 shows more detail of S302 shown in Fig. 19. 

Fig. 22 shows more detail of S308 shown in Fig. 19. 

Fig. 23 shows an image based on data for a tornado at 
time t=tl. 

Fig. 24 shows an image based on data for a tornado at 
10 time t=t2. 

Fig. 25 shows an intermediate image obtained based on 



ill 

py the images of Figs. 23 and 24. 

h 
m 

^ 15 DETAILED DESCRIPTION OF THE INVENTION 



The invention will now be described based on the 
preferred embodiments, which are not intended to limit the 
scope but to exemplify the present invention. All of the 
.20 features and the combinations thereof described in the 

embodiments are not necessarily essential to the invention. 

First, the multiresolutional critical point filter 
technology and the image matching processing using the 
technology, both of which will be utilized in the preferred 
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embodiments, will be described in detail as "Base Technology''. 
Namely, the following sections [1] and [2] (below) belong to 
the base technology, where section [1] describes elemental 
techniques and section [2] describes a processing procedure. 
5 These techniques are patented under Japanese Patent No. 
2 927350 and owned by the same assignee of the present 
invention. As described in more detail below following the 
discussion of the base technology, according to the 
embodiments of the present invention there is provided a mesh 
;;J 10 on an image, so that lattice points of the mesh represent a 
i| plurality of pixels of the image. Thus, even though 

ry application efficiency for a pixel-by-pixel matching technique 

O described in the base technology is naturally high, it is 

iU 

to be noted that the image matching techniques provided in the 
!:| 15 present embodiments are not limited to the same levels. 

After discussing the base technology, there is a 
description, related to Figs. 18 to 23, of image data coding 
and decoding techniques for multivariate space processing 
according to embodiments of the present invention and 
20 utilizing the base technology. 

Base Technology 

[1] Detailed description of elemental techniques 
[1.1] Introduction 
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Using a set of new multiresolutional filters called 
critical point filters, image matching is accurately computed. 
There is no need for any prior knowledge concerning the 
content of the images or objects in question. The matching of 
5 the images is computed at each resolution while proceeding 
through the resolution hierarchy. The resolution hierarchy 
proceeds from a coarse level to a fine level. Parameters 
necessary for the computation are set completely automatically 

M by dynamical computation analogous to human visual systems. 

;<3 10 Thus, There is no need to manually specify the correspondence 
of points between the images. 

i 

ru '^h® hase technology can be applied to, for instance, 

completely automated morphing, object recognition, stereo 

ru 

photogrammetry, volume rendering, and smooth generation of 
Q 15 motion images from a small number of frames. When applied to 
morphing, given images can be automatically transformed. When 
applied to volume rendering, intermediate images between cross 
sections can be accurately reconstructed, even when a distance 
between cross sections is rather large and the cross sections 
20 vary widely in shape. 

[1.2] The hierarchy of the critical point filters 

The multiresolutional filters according to the base 
technology preserve the intensity and location of each 
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critical point included in the images while reducing the 
resolution. Initially, let the width of an image to be 
examined be N and the height of the image be M. For 
simplicity, assume that N=M=2n where n is a positive integer. 
An interval [0, N] c R is denoted by I . A pixel of the image 
at position {i, j) is denoted by p<^'^' where i,j e I. 

Here, a multiresolutional hierarchy is introduced. 
Hierarchized image groups are produced by a multiresolutional 
filter. The multiresolutional filter carries out a two 
dimensional search on an original image and detects critical 
points therefrom. The multiresolutinal filter then extracts 
the critical points from the original image to construct 
another image having a lower resolution. Here, the size of 
each of the respective images of the m-th level is denoted as 
2"X2"' (0<m<n) . A critical point filter constructs the 
following four new hierarchical images recursively, in the 
direction descending from n. 

— (1) 

where we let 

P(,ij) -Poj) -P(U) -P(.U) -P(ij) <2) 
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The above four images are referred to as subimages 
hereinafter. When minx^t^x+i and maxx^ts:x+i are abbreviated to a 
and P, respectively, the subimages can be expressed as 
follows : 

Namely, they can be considered analogous to the tensor 
products of a and p. The subimages correspond to the 
respective critical points. As is apparent from the above 
equations, the critical point filter detects a critical point 
of the original image for every block consisting of 2 X 2 
pixels. In this detection, a point having a maximum pixel 
value and a point having a minimum pixel value are searched 
with respect to two directions, namely, vertical and 
horizontal directions, in each block. Although pixel 
intensity is used as a pixel value in this base technology, 
various other values relating to the image may be used. A 
pixel having the maximum pixel values for the two directions, 
one having minimum pixel values for the two directions, and 
one having a minimum pixel value for one direction and a 
maximum pixel value for the other direction are detected as a 
local maximum point, a local minimum point, and a saddle 
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point, respectively. 

By using the critical point filter, an image (1 pixel 
here) of a critical point detected inside each of the 
respective blocks serves to represent its block image (4 
5 pixels here) in the next lower resolution level. Thus, the 
resolution of the image is reduced. From a singularity 
theoretical point of view, a{x)a(y) preserves the local 
minimiom point (minima point ), p (x) P (y) preserves the local 
maximum point (maxima point) , a(x)p(y) andp(x)a(y) preserve 

3 

\| 10 the saddle points. 

m 

•M At the beginning, a critical point filtering process is 

iU applied separately to a source image and a destination image 

n 

III which are to be matching-computed. Thus, a series of image 

J groups, namely, source hierarchical images and destination 

py 15 hierarchical images aire generated. Four source hierarchical 

images and four destination hierarchical images are generated 

corresponding to the types of the critical points. 

Thereafter, the source hierarchical images and the 

destination hierarchical images are matched in a series of 
20 resolution levels. First, the minima points are matched using 

p(m,o)_ Next, the first saddle points are matched using p'"''^' 

based on the previous matching result for the minima points. 

The second saddle points are matched using p*""'^'. Finally, 

the maxima points are matched using p^""'^'. 



Figs. Ic and Id show the subimages p<5'°' of the images 
in Figs, la and lb, respectively. Similarly, Figs, le and If 
show the subimages p'^'^'. Figs. Ig and Ih show the subimages 
p^^'^^ and Figs, li and Ij show the subimages p<^'^'. 
Characteristic parts in the images can be easily matched using 
subimages. The eyes can be matched by p'^'°> since the eyes 
are the minima points of pixel intensity in a face. The 
mouths can be matched by p<^'^' since the mouths have low 
intensity in the horizontal direction. Vertical lines on both 
sides of the necks become clear by p'^'^'. The ears and bright 
parts of the cheeks become clear by p'^'^' since these are the 
maxima points of pixel intensity. 

As described above, the characteristics of an image can 
be extracted by the critical point filter. Thus, by 
comparing, for example, the characteristics of an image shot 
by a camera with the characteristics of several objects 
recorded in advance, an object shot by the camera can be 
identified. 

[1.3] Computation of mapping between images 

Now, for matching images, a pixel of the source image at 
the location (i,j) is denoted by /7f"J.^ and that of the 
destination image at (k,l) is denoted by ^{g^ where ,i, j, k, 1 
e I. The energy of the mapping between the images (described 
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later in more detail) is then defined. This energy is 
determined by the difference in the intensity of the pixel of 
the source image and its corresponding pixel of the 
destination image and the smoothness of the mapping. First, 
5 the mapping f (-'O :p(-.o) ^ g(ra,o) between p''"'"' and q'^-'^' with the 
minimum energy is computed. Based on f^'^'°\ the mapping f(™'i> 
between p'^'^' and q^™'^' with the minimum energy is computed. 
This process continues until f'"''^) between p<'*'35 and qf'"'^) j_g 
computed. Each f'""'^' (i = 0,1,2,...) is referred to as a 
10 submapping. The order of i will be rearranged as shown in the 
following equation (3) in computing f<™'^' for reasons to be 
described later. 

y (»»>') . pim.trO)) ^ qimMl)) P J 

where a(i) e {0,1,2,3}. 

15 

[1. 3. 1] Bijectivity 

When the matching between a source image and a 
destination image is expressed by means of a mapping, that 
mapping shall satisfy the Bijectivity Conditions (BC) between 
20 the two images (note that a one-to-one surjective mapping is 
called a bijection) . This is because the respective images 
should be connected satisfying both surjection and injection, 
and there is no conceptual supremacy existing between these 
images. It is to be noted that the mappings to be constructed 
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here are the digital version of the bijection. In the base 
technology, a pixel is specified by a co-ordinate point. 

The mapping of the source subimage (a subimage of a 
source image) to the destination subimage (a subimage of a 
destination image) is represented by f '""'S' : 1/2''""' X 1/2"""' 
1/2"-" X 1/2"-=" (s = 0,1,...), where /^gf means that pf^g of 

.the source image is mapped to q'^^f^ of the destination image. 

For simplicity, when f (i,j) = {k,l) holds, a pixel q(k,i) is 

denoted by qf(i,j). 

When the data sets are discrete as image pixels (grid 

points) treated in the base technology, the definition of 

bijectivity is important. Here, the bijection will be defined 

in the following manner, where i, j, k and 1 are all integers. 

First, a square region R defined on the source image plane is 

considered 

P^ij) Pii+ij}Po+i,j+i)P(u+i) 

where i = 0, 2"-!, and j = 0, 2^-1. The edges of R are 
directed as follows: 

P(U) P0*UJ)'P(MJ)P(MJ+1)'P(M,M)P(U*1) P0,M)P(U) (5) 

This square region R will be mapped by f to a 
quadrilateral on the destination image plane: 

^f(U)^%Xj)^%+ij+i)'i%j+i) (6) 

This mapping f <"''^' (R) , that is, 
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should satisfy the following bijectivity conditions (referred 
to as BC hereinafter) : 

1. The edges of the quadrilateral f ^""'^^ (R) should not 
intersect one another. 

2. The orientation of the edges of f'^^'^MR) should be the same 
as that of R (clockwise in the case shown in Fig. 2, described 
below) . 

3. As a relaxed condition, a retraction mapping is allowed. 

Without a certain type of a relaxed condition as in, for 
example, condition 3 above, there would be no mappings which 
completely satisfy the BC other than a trivial identity 
mapping. Here, the length of a single edge of f (r) ^ay be 
zero. Namely, f '"''^^ (R) may be a triangle. However, f (R) 
is not allowed to be a point or a line segment having area 
zero. Specifically speaking, if Fig. 2R is the original 
quadrilateral. Figs. 2A and 2D satisfy the BC while Figs 2B, 
2C and 2E do not satisfy the BC. 

In actual implementation, the following condition may be 
further imposed to easily guarantee that the mapping is 
surjective. Namely, each pixel on the boundary of the source 
image is mapped to the pixel that occupies the same location 
at the destination image. In other words, f(i,j)=(i,j) (on 
the four lines of i=0, i=2"^-l, j=0, j=2"^-l) . This condition 
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will be hereinafter referred to as an additional condition. 



[1. 3. 2] Energy of mapping 

[1. 3. 2. 1] Cost related to the pixel intensity 
5 The energy of the mapping f is defined. An objective 

here is to search a mapping whose energy becomes minimum. The 
energy is determined mainly by the difference in the intensity 
between the pixel of the source image and its corresponding 
Q pixel of the destination image. Namely, the energy Cj"j*^ of 

%^ 10 the mapping f'"''^' at(i,j) is determined by the following 
equation (7) . 

;^ ^£?=h/'fD?)-f^(^)^3))r — (7) 

I,': where ) and are the intensity values of the 

II pixels pl^jf and q%fj^, respectively. The total energy 0^™'=^ of 

15 f is a matching evaluation equation, and can be defined as the 
sum of C^^jf as shown in the following equation (8) . 

/=0 7=0 

[1. 3. 2. 2] Cost related to the locations of the pixel for 
20 smooth mapping 

In order to obtain smooth mappings, another energy Df 
for the mapping is introduced. The energy Df is determined by 
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the locations of p'^J'j'^ and g-J^/J^ ( i=0 , 1 , 2"'-l , j=0/l, 2"-!) , 
regardless of the intensity of the pixels. The energy D^^f^ of 
the mapping f'""'^' at a point (i,j) is determined by the 
following equation (9) . 

where the coefficient parameter r\ which is equal to or greater 
than 0 is a real number. And we have 

Eo^^=\\(U)-P"'''\U)f — (10) 

= t II 1^4 — (11) 

where 

\\(x,y)\\ = ^x'+y' (12), 

i' and j' are integers and f(i',j') is defined to be zero for 
i'<0 and j''<0. Eo is determined by the distance between (i,j) 
and f(i,j), Eo prevents a pixel from being mapped to a pixel 
too far away from it. However, as explained below, Eo can be 
replaced by another energy function. Ei ensures the 
smoothness of the mapping. Ei represents a distance between 
the displacement of p(i,j) and the displacement of its 
neighboring points. Based on the above' consideration, another 
evaluation equation for evaluating the matching, or the energy 
Df is determined by the following equation: 
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/=2''-lj=2"'-l 

/)<"-)= 2 Z^&," — (13) 

[1. 3. 2. 3] Total energy of the mapping 

The total energy of the mapping, that is, a combined 
5 evaluation equation which relates to the combination of a 

plurality of evaluations, is defined as Klf-"^ +Df^^ , where 9^0 
is a real number. The goal is to detect a state in which the 
p combined evaluation equation has an extreme value, namely, to 

find a mapping which gives the minimum energy expressed by the 

m 

»R 10 following: 

min{AC}»''*>+Z)f'^>} (14) 

lb 

it-l Care must be exercised in that the mapping becomes an 

identity mapping if >.=0 and r|=0 (i.e., f (i, j ) = (i, j ) for all 

^'■^ i=0,l,...,2'"-l and j=0 , 1 , 2'"-l) . As will be described later, 

15 the mapping can be gradually modified or transformed from an 
identity mapping since the case of X=0 and ti=0 is evaluated at 
the outset in the base technology. If the combined evaluation 
equation is defined as Cf'^'^^mf''^ where the original position 
of X is changed as such, the equation with 1=0 and r[=0 will be 
20 Cj"'*^ only. As a result thereof, pixels would randomly matched 
to each other only because their pixel intensities are close, 
thus making the mapping totally meaningless. Transforming the 
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mapping based on such a meaningless mapping makes no sense. 
Thus, the coefficient parameter is so determined that the 
identity mapping is initially selected for the evaluation as 
the best mapping. 

Similar to this base technology, differences in the 
pixel intensity and smoothness are considered in a technique 
called optical flow" that is known in the art. However, the 
optical flow technique cannot be used for image transformation 
since the optical flow technique takes into account only the 
local movement of an object. However, global correspondence 
can also be detected by utilizing the critical point filter 
according to the base technology. 

[1. 3. 3] Determining the mapping with mult iresolut ion 
A mapping fnu.n which gives the minimiam energy and 
satisfies the BC is searched by using the multiresolution 
hierarchy. The mapping between the source subimage and the 
destination subimage at each level of the resolution is 
computed. Starting from the top of the resolution hierarchy 

(i.e., the coarsest level), the mapping is determined at each 
resolution level, and where possible, mappings at other levels 
are considered. The number of candidate mappings at each 
level is restricted by using the mappings at an upper (i.e., 
coarser) level of the hierarchy. More specifically speaking, 
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in the course of determining a mapping at a certain level, the 
mapping obtained at the coarser level by one is imposed as a 
sort of constraint condition. 

We thus define a parent and child relationship between 
resolution levels. When the following equation (15) holds, 

where [xj denotes the largest integer not exceeding x, 
PirJ'f and g^^jl;'^ are respectively called the parents of p^^^f 
and ^Ti?'- Conversely, p^^^^f and q^^j^^ are the child of p^J^^ 
and the child of gl^jY\ respectively. A function parent (i,j) 
is defined by the following equation (16) : 



parent(i,j) = (\ 




— (16) 



Now, a mapping between /^f")"^ and q^^jf is determined by 
computing the energy and finding the minimum thereof. The 
value of f '""'^^ (i, j) = (k,l) is determined as follows using f (m- 
l,s) (m=l,2,...,n) . First of all, a condition is imposed that 
^(kf) should lie inside a quadrilateral defined by the 
following definitions (17) and (18). Then, the applicable 
mappings are narrowed down by selecting ones that are thought 
to be reasonable or natural among them satisfying the BC. 
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where 



g^'"'"\i,j) = /^""^-"Kparentii,])) + f^"'-'''\parent(i,j) + (1,1)) — { 18 ) 
The quadrilateral defined above is hereinafter referred 
to as the inherited quadrilateral of p^^'^ . The pixel 
minimizing the energy is sought and obtained inside the 
inherited quadrilateral. 

Fig. 3 illustrates the above-described procedures. The 
pixels A, C and D of the source image are mapped to A' , B' , 
C and D' of the destination image, respectively, at the (m- 
l)th level in the hierarchy. The pixel pj";^^ should be mapped 
to the pixel ^^wj,^^) which exists inside the inherited 
quadrilateral A'B'C'D'. Thereby, bridging from the mapping at 
the (m-l)th level to the mapping at the m-th level is 
achieved. 

The energy Eq defined above may now be replaced by the 
following equations (19) and (20) : 



for computing the submapping f*"''^' and the submapping f*"''^' at 
the m-th level, respectively. 

In this manner, a mapping which maintains a low energy 
of all the submappings is obtained. Using the equation (20) 




— (19) 



= Wf-'-'a, J) - f""''-'\i, j)f, (1 < 0 — (20) 



MN-70022 



makes the submappings corresponding to the different critical 
points associated to each other within the same level in order 
that the subimages can have high similarity. The equation 
(19) represents the distance between f""'^'(i,j) and the 
location where (i,j) should be mapped when regarded as a part 
of a pixel at the (m-l)the level. 

When there is no pixel satisfying the BC inside the 
inherited quadrilateral A' B' C D' , the following steps are 
taken. First, pixels whose distance from the boundary of 
A'D'C'D' is L (at first, L=l) are examined. If a pixel whose 
energy is the minimum among them satisfies the BC, then this 
pixel will be selected as a value of f ^'^'^^ (i, j ) . L is 
increased until such a pixel is found or L reaches its upper 
bound L^^l^. L'-^^ is fixed for each level m. If no pixel is 
found at all, the third condition of the BC is ignored 
temporarily and such mappings that caused the area of the 
transformed quadrilateral to become zero (a point or a line) 
will be permitted so as to determine f*"'^'(i,j). If such a 
pixel is still not found, then the first and the second 
conditions of the BC will be removed. 

Multiresolution approximation is essential to 
determining the global correspondence of the images while 
preventing the mapping from being affected by small details of 
the images. Without the multiresolution approximation, it is 
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impossible to detect a correspondence between pixels whose 
distances are large. In the case where the multiresolution 
approximation is not available, the size of an image will 
generally be limited to a very small size, and only tiny 
5 changes in the images can be handled. Moreover, imposing 

smoothness on the mapping usually makes it difficult to find 
the correspondence of such pixels. That is because the energy 
of the mapping from one pixel to another pixel which is far 
therefrom is high. On the other hand, the multiresolution 
10 approximation enables finding the approximate correspondence 
of such pixels- This is because the distance between the 
pixels is small at the upper (coarser) level of the hierarchy 
of the resolution. 

15 [1. 4] Automatic determination of the optimal parameter values 
One of the main deficiencies of the existing image 
matching techniques lies in the difficulty of parameter 
adjustment. In most cases, the parameter adjustment is 
performed manually and it is extremely difficult to select the 

20 optimal value. However, according to the base technology, the 
optimal parameter values can be obtained completely 
automatically. 

The systems according to this base technology include 
two parameters, namely, X and ti, where X, and r\ represent the 



weight of the difference of the pixel intensity and the 
stiffness of the mapping, respectively. In order to 
automatically determine these parameters, the are initia-lly 
set to 0. First, % is gradually increased from X=0 while r\ is 
fixed at 0. As X becomes larger and the value of the combined 
evaluation equation (equation (14)) is minimized, the value of 
Cy^*^ for each submapping generally becomes smaller. This 
basically means that the two images are matched better. 
However, if X exceeds the optimal value, the following 
phenomena occur: 

1. Pixels which should not be corresponded are erroneously 
corresponded only because their intensities are close. 

2. As a result, correspondence between images becomes 
inaccurate, and the mapping becomes invalid. 

3. As a result, Dj""'^ in equation (14) tends to increase 
abruptly. 

4. As a result, since the value of equation (14) tends to 
increase abruptly, f*"*'^' changes in order to suppress the 
abrupt increase of D)"'*^ . As a result, C}"'*^ increases. 

Therefore, a threshold value at which C}"'*' turns to an 
increase from a decrease is detected while a state in which 
equation (14) takes the minimum value with X being increased 
is kept. Such X is determined as the optimal value at ti=0 . 

MN-70022 



32 

Next, the behavior of Cf'"'* is examined while r| is increased 
gradually, and ri will be automatically determined by a method 
described later, ' X will then again be determined 
corresponding to such an automatically determined ti. 

The above-described method resembles the focusing 
mechanism of human visual systems. In the human visual 
systems, the images of the respective right eye and left eye 
are matched while moving one eye. When the objects are 
clearly recognized, the moving eye is fixed. 

[1. 4. 1] Dynamic determination o^.X 

Initially, X is increased from 0 at a certain interval, 
and a subimage is evaluated each time the value of X changes. 
As shown in equation (14), the total energy is defined by 
^^(».*) _^ ^ equation (9) represents the smoothness 

and theoretically becomes minimum when it is the identity 
mapping. Eq and Ei increase as the mapping is further 
distorted. Since Ei is an integer, 1 is the smallest step of 
Df-'K Thus, it is impossible to change the mapping to reduce 
the total energy unless a changed amount (reduction amount) of 
the current ACf/"j)^ is equal to or greater than 1. Since Df'^^ 
increases by more than 1 accompanied by the change of the 
.mapping, the total energy is not reduced unless ACjf))^ is 
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reduced by more than 1. 

Under this condition, it is shown that Cjl^jf decreases in 
normal cases as X increases. The histogram of q'Jf is denoted 
as h(l), where h(l) is the number of pixels whose energy C^^j^^ 
is 1^. In order that Xl^ ^ 1, for example, the case of 1^=1 A 
is considered. When X varies from Xi to ^2, a number of pixels 
(denoted A) expressed by the following equation (21) : 



i = y /.(/) ^ / hwi = -]h(i)^A = )f 

Kil '-i ^ ^'^ 



(21) 



changes to a more stable state having the energy shown in 
equation (22) : 

Cr^-/^=Cr-^-}. — (22) 

Here, it is assumed that the energy of these pixels is 



changes by: 

5Cr>=-| — (23) 

As a result, equation (24) holds. 

dC^/-'^ h(l) 

= — (24) 

Since h(l)>0 , C)'"'*^ decreases in the normal case. However, 
when X exceeds the optimal value, the above phenomenon, that 



is, an increase in Cf'"'^ occurs. The optimal value of X, is 
determined by detecting this phenomenon. 
When 

hQ) = Hl'=^ (25) 

is assumed, where both H(H>0) and k are constants, the 
equation (26) holds: 

dA Z''^^"^ ^^^^ 
Then, if k^t-S, the following equation (27) holds: 

C/'"'=C+ (3/2 + ^^)^3™ — (27) 

The equation (27) is a general equation of C|'"'*^ .(where C is a 
constant) . 

When detecting the optimal value of X, the number of 
■ pixels violating the BC may be examined for safety. In the 
course of determining a mapping for each pixel, the 
probability of violating the BC is assumed as a value po here. 
In this case, since 
dA h(l) 

holds, the nijmber of pixels violating the BC increases at a 
rate of: 

_ h(l)p^ 
^0 — Tiir- — (29) 



Thus, 

-^ = 1 (30) 

is a constant. If it is assumed that h(l)=Hl'', the following 
equation (31), for example, 

B,X"'^"'=p,H (31) 

becomes a constant. However, when X exceeds the optimal 
value, the above value of equation (31) increases abruptly. 
By detecting this phenomenon, i.e. whether or not the value of 
BqJ^''^*'''^ 12'" exceeds an abnormal value B^^^^^^ , the optimal value 
of X can be determined. Similarly, whether or not the value 
of 5,.i^'^**'V2'" exceeds an abnormal value B^^^^^ can be used to 
check for an increasing rate Bi of pixels violating the third 
condition of the BC. The reason why the factor 2"" is 
introduced here will be described at a later stage. This 
system is not sensitive to the two threshold values 5^^^^^^ and 
The two threshold values ^o,^^^ and B^^^^^ can be used to 
detect excessive distortion of the mapping which may not be 
detected through observation of the energy C^'"'*> . 

In the experimentation, when X exceeded 0.1 the 
computation of f = ' was stopped and the computation of f<i"'S+i) 
was started. That is because the computation of submappings 
is affected by a difference of only 3 out of 255 levels in 
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pixel intensity when X > 0.1 and it is then difficult to 
obtain a correct result. 

[1. 4. 2] Histogram h(l) 

The examination of C^/''^ does not depend on the histogram 
h(l), however, the examination of the BC and its third 
condition may be affected by h(l). When {X , Cf''^) is actually 
plotted, k is usually close to 1. In the experiment, k=l is 
used, that is, 8^2^ and B^^^ are examined. If the true value 
of k is less than 1, B^^ and B^^ are not constants and 
increase gradually by a factor of A^'"*^'^ . If h(l) is a 
constant, the factor is, for example, ^'^ . However, such a 
difference can be absorbed by setting the threshold B^^^ 
appropriately. 

Let us model the source image by a circular object, with 
its center at(xo,yo) and its radius r, given by: 



P(ij) 



0... (otherwise) 



— (32) 



and the destination image given by: 
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— (33) 

with its center at(xi,yi) and radius r. In the above, let 
c{x) have the form of c(x)=x''. When the centers {xo,yo) and 
(xi,yi) are sufficiently far from each other, the histogram 
h(l) is then in the form: 

hiD'Krl'ik^Q) (34) 

When k=l, the images represent objects with clear 
boundaries embedded in the background. These objects become 
darker toward their centers and brighter toward their 
boundaries. When k=-l, the images represent objects with 
vague boundaries. These objects are brightest at their 
centers, and become darker toward their boundaries. Without 
much loss of generality, it suffices to state that objects in 
images are generally between these two types of objects. 
Thus, choosing k such that -l<k<l can cover most cases and the 
equation (27) is generally a decreasing function for this 
range . 

As can be observed from the above equation (34), 
attention must be directed to the fact that r is influenced by 
the resolution of the image, that is, r is proportional to a"'. 
This is the reason for the factor 2"^ being introduced in the 
above section [1.4.1]. 

[1. 4. 3] Dynamic determination of ti 
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The parameter T| can also be automatically determined in 
a similar manner. Initially, ti is set to zero, and the final 
mapping f''^' and the energy Cf^ at the finest resolution are 
computed. Then, after ri is increased by a certain value A77 , 
the final mapping f^''' and the energy C}"^ at the finest 
resolution are again computed. This process is repeated until 
the optimal value of t| is obtained. x\ represents the 
stiffness of the mapping because it is a weight of the 
following equation (35) : 

Eot^ =\f""''\iJ)-f""'-'\Uj)f — (35) 

If T| is zero, Df^ is determined irrespective of the 
previous submapping, and the present submapping may be 
elastically deformed and become too distorted. On the other 
hand, if t] is a very large value, Df is almost completely 
determined by the immediately previous submapping. The 
submappings are then very stiff, and the pixels are mapped to 
almost the same locations. The resulting mapping is therefore 
the identity mapping. When the value of t] increases from 0, 
Cj"^ gradually decreases as will be described later. However, 
when the value of rj exceeds the optimal value, the energy 
starts increasing as shown in Fig. 4. In Fig. 4, the x-axis 
represents r\, and y-axis represents Cf. 
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The optimum value of if] which minimizes C^^ can be 
obtained in this manner. However, since various elements 
affect this computation as compared to the case of 1, Cf^ 
changes while slightly fluctuating. This difference is caused 
because a subraapping is re-computed once in the case of X, 
whenever an input changes slightly, whereas all the 
submappings must be re-computed in the case of t|. Thus, 
whether the obtained value of Cf^ is the minimum or not cannot 
be determined as easily. When candidates for the minimum 
value are found, the true minimum needs to be searched by 
setting up further finer intervals. 

[1. 5] Super sampling 

When deciding the correspondence between the pixels, the 
range of f can be expanded to R X R (R being the set of 
real numbers) in order to increase the degree of freedom. In 
this case, the intensity of the pixels of the destination 
image is interpolated, to provide f'""'^' having an intensity at 
non-integer points: 

»'(*?4„^) — (36) 

That is, supersampling is performed. In an example 
implementation, = ) may take integer and half integer 
values, and 
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^(^S)W5)) (37) 

is given by 

(^(9f;::;?)+^(^(?);!<M)))/2 — os) 

5 [ 1. 6] Normalization of the pixel intensity of each image 
When the source and destination images contain quite 
different objects, the raw pixel intensity may not be used to 
compute the mapping because a large difference in the pixel 
q intensity causes excessively large energy C^J"'"^ and thus making 

Cn 10 it difficult to obtain an accurate evaluation. 

For example, a matching between a human face and a cat's 
face is computed as shown in Figs. 20(a) and 20(b). The cat's 
\^ face is covered with hair and is a mixture of very bright 

Q pixels and very dark pixels. In this case, in order to 

m 

15 compute the submappings of the two faces, subimages are 

normalized. That is, the darkest pixel intensity is set to 0 
while the brightest pixel intensity is set to 255, and other 
pixel intensity values are obtained using linear 
interpolation . 

20 

[ 1 . 7 ] Implementation 

In an example implementation, a heuristic method is 
utilized wherein the computation proceeds linearly as the 
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source image is scanned. First, the value of f is 
determined at the top leftmost pixel (i,j)=(0,0). The value 
of . each f^"''^'(i,j) is then determined while i is increased by 
one at each step. When i reaches the width of the image, j is 
increased by one and i is reset to zero. Thereafter, 
f^°''^^(i,j) is determined while scanning the source image. 
Once pixel correspondence is determined for all the points, it 
means that a single mapping f''^'^^ is determined. 

When a corresponding point qf(i,j) is determined for P(i,j), 
a corresponding point qf(i,j+i) of P(i,j+i) is determined next. 
The position of qf(i,j+i) is constrained by the position of qf(i,j) 
since the position of qf(x,j+i) satisfies the BC. Thus, in this 
system, a point whose corresponding point is determined 
earlier is given higher priority. If the situation continues 
in which (0,0) is always given the highest priority, the final 
mapping might be unnecessarily biased. In order to avoid this 
bias, f'"''^' is determined in the following manner in the base 
technology. 

First, when (s mod 4) is 0, f'""'^' is determined starting 
from (0,0) while gradually increasing both i and j. When (s 
mod 4) is 1, f*"''^' is determined starting from the top 
rightmost location while decreasing i and increasing j . When 
(s mod 4) is 2, f'""'^' is determined starting from the bottom 
rightmost location while decreasing both i and j. When (s mod 
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4) is 3, f^'"'^' is determined starting from the bottom leftmost 
location while increasing i and decreasing j . Since a concept 
such as the submapping, that is, a parameter s, does not exist 
in the finest n-th level, f*"*'^^ is computed continuously in two 
directions on the assumption that s=0 and s=2 . 

In this implementation, the values of f*"''^'(i,j) 
(m=0,...,n) that satisfy the BC are chosen as much as possible 
from the candidates (k,l) by imposing a penalty on the 
candidates violating the BC. The energy D(k,i) of a candidate 
that violates the third condition of the BC is multiplied by (p 
and that of a candidate that violates the "first or second 
condition of the BC is multiplied by In this 

implementation, (p=2 and \|/=100000 are used. 

In order to check the above-mentioned BC, the following 
test may be performed as the procedure when determining 
(k,l)=f ^'^'^^ (i, j) . Namely, for each grid point (k, 1) in the 
inherited quadrilateral of f (i, j ) , whether or not the z- 
component of the outer product of 

W = AxB (39) 

is equal to or greater than 0 is examined, where 

^ = *<5t!,„_„?<Z>.<,,,.„ —(40, 

5=«<Z!,„.„«S;? —(41) 

Here, the vectors are regarded as 3D vectors and the z-axis is 
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defined in the orthogonal right-hand coordinate system. When 
W is negative, the candidate is imposed with a penalty by 
multiplying Dg;*> by \|/ so that it is not as likely to be 
selected. 

5 Figs. 5(a) and 5(b) illustrate the reason why this 

condition is inspected. Fig. 5(a) shows a candidate without a 
penalty and Fig. 5(b) shows one with a penalty. When 
determining the mapping (i, j+1) for the adjacent pixel at 

■Q (i,j+l), there is no pixel on the source image plane that 

S 10 satisfies the BC if the z-component of W is negative because 

m 

•\P then q^^f^^ passes the boundary of the adjacent quadrilateral. 

h 

rU [1. 7. 1] The order of submappings 

In this implementation, ct(0)=0, a(l)=l, o(2)=2, o(3)=3, 
15 a(4)=0 are used when the resolution level is even, while 
a(0)=3, <j(l)=2, a(2)=l, a(3)=0, a(4)=3 are used when the 
resolution level is odd. Thus, the submappings are shuffled 
to some extent. It is to be noted that the submappings are 
primarily of four types, and s may be any of 0 to 3. However, 
20 a processing with s=4 is used in this implementation for a 
reason to be described later. 



[1. 8] Interpolations 

After the mapping between the source and destination 
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images is determined, the intensity values of the 

corresponding pixels are interpolated. In the implementation, 

trilinear interpolation is used. Suppose that a square 

P<i,j)P(i+i,j)P(i+i, j+i)P(i,j+i) on the source image plane is mapped to 

a quadrilateral qf (i, j)qf (i+i, j)qf (i+i, j+i)qf (i, j+i) on the destination 

image plane. For simplicity, the distance between the image 

planes is assumed to be 1. The intermediate image pixels 

r(x,y,t) (0^x<N-l, 0:^y^M-l) whose distance from the source 

image plane is t (0<t<l) are obtained as follows. First, the 

location of the pixel r(x,y,t), where x,y,tGR, is determined 

by equation (42) : 

(x, y) = il- dxXl - dyXl - t)ii, J) + (1 - dxXl - dy)tf{i, j) 
+ dxil - dyX\ - t)(i + 1, j) + dxO- - dy)tf(i + 1, J) 
+ (l-dx)dy(l-t)(i,j + l) + (l-dx)dytf(i,J + l) ~~~ ^^^^ 
+ dxdy(l - t)(i + 1, y + 1) + dxdytf(i + l,j+ 1) 

The value of the pixel intensity at r(x,y,t) is then. 

determined by equation (43) : i 

V(r(x, y, 0) = (1 - dx)(l - dy)(l - t)V(p^^ .^ ) + (1 - ch)(l - dy)tV(q^^, j^) 



where dx and dy are parameters varying from 0 to 1 . 

[1. 9] Mapping to which constraints are imposed 

So far, the determination of a mapping in which no 
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+ dx(l - dyXl - tWip^urj^ ) + dx(l- dy)tV{q^ 
+ (1 - dx)dy{l - tW(Pcu-,y ) + (1 - dx)dytV(q^ 

+ dxdy(l - 0J^(;?(,-., ,y.,) ) + ) 




— (43) 



constraints are imposed has been described. However, if a 
correspondence between particular pixels of the source and 
destination images is provided in a predetermined manner, the 
mapping can be determined using such correspondence as a 
constraint . 

The basic idea is that the source image is roughly 
deformed by an approximate mapping which maps the specified 
pixels of the source image to the specified pixels of the 
destination image and thereafter a mapping f is accurately 
computed. 

First, the specified pixels of the source image are 
mapped to the specified pixels of the destination image, then 
the approximate mapping that maps other pixels of the source 
image to appropriate locations are determined. In other words, 
the mapping is such that pixels in the vicinity of a specified 
pixel are mapped to locations near the position to which the 
specified one is mapped. Here, the approximate mapping at the 
m-th level in the resolution hierarchy is denoted by F*""^ . 

The approximate mapping F is determined in the following 
manner. First, the mappings for several pixels are specified. 
When ns pixels 

PQo,JoXP(h,j\\-,P(L^-iJr.,,-i) — (44) 

of the source image are specified, the following values in the 
equation (45) are determined. 
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F"'\h,j\) = (f^rJr),-, (45) 

For the remaining pixels of the source image, the amount 
of displacement is the weighted average of the displacement of 
P(ih/jh) (h=0,..., ng-l) . Namely, a pixel P(i,j) is mapped to the 
following pixel (expressed by the equation (46)) of the 
destination image. 

(hJ)+ -h,h -jH>eight^ii,j) 
F'-\iJ) = '-^ — (46) 



total _ weight{i, j) 



total _weight(i,j) = 'X'l/IK/, -jf 



(48 



Second, the energy D^^j^^ of the candidate mapping f is 
changed so that a mapping f similar to f'""' has a lower energy. 

where 



I '-J 

[IfC") (/, J) - /("-^^ (/, ;)f , otherwise 



(50) 



where K,p>0. Finally, the resulting mapping f is determined 
by the above-described automatic computing process. 

Note that E^'^-^J becomes 0 if f<"''^'(i,j) is sufficiently 
close to F<""(i,j) i.e., the distance therebetween is equal to 
or less than 




— (51) 



This has been defined in this way because it is desirable to 
determine each value f'"''^'(i,j) automatically to fit in an 
appropriate place in the destination image as long as each 
value f^'"'^'(i,j) is close to F'""(i,j), For this reason, there 
is no need to specify the precise correspondence in detail to 
have the source image automatically mapped so that the source 
image matches the destination image. 

[2] Concrete Processing Procedure 

The flow of a process utilizing the respective elemental 
techniques described in [1] will now be described. 

Fig. 6 is a flowchart of the overall procedure of the 
base technology. Referring to Fig. 6, a source image and 
destination image are first processed using a 

multiresolutional critical point filter (SI) . The source image 
and the destination image are then matched (S2) . As will be 
understood, the matching (S2) is not required in every case, 
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and other processing such as image recognition may be 
performed instead, based on the characteristics of the source 
image obtained at SI. 

Fig. 7 is a flowchart showing details of the process SI 
shown in Fig. 6. This process is performed on the assumption 
that a source image and a destination image are matched at S2. 
Thus, a source image is first hierarchized using a critical 
point filter (SIO) so as to obtain a series of source 
hierarchical images. Then, a destination image is hierarchized 
in the similar manner (SIX) so as to obtain a series of 
destination hierarchical images. The order of SIO and Sll in 
the flow is arbitrary, and the source image and the 
destination image can be generated in parallel. It may also be 
possible to process a number of source and destination images 
as required by subsequent processes. 

Fig. 8 is a flowchart showing details of the process at 
SIO shown in Fig. 7. Suppose that the size of the original 
source image is 2"X2". Since source hierarchical images are 
sequentially generated 'from an image with a finer resolution 
to one with a coarser resolution, the parameter m which 
indicates the level of resolution to be processed is set to n 
(SlOO) . Then, critical points are detected from the images 
P'™'°^' P^"'''\ p'"'^' and p<"'^>of the m-th level of resolution, 
using a critical point filter (SlOl) , so that the images p^"*" 
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1,0,^ p(m-i,i)^ p(m-i,2, p(n.-i,3) (ra-l)th level are 

generated (S102) . Since m=n here, p<™'°' =p("''^> =p<"''2) ^p(ni,3) 
=p'''' holds and four types of sub images are thus generated from 
a single source image. 
5 Fig. 9 shows correspondence between partial images of 

the m-th and those of (m-l)th levels of resolution. Referring 
to Fig. 9, respective numberic values shown in the figure 
represent the intensity of respective pixels. p'^^-^' 
|;J symbolizes any one of four images p<"'°' through p'"''^', and when 

10 generating p<ia-i,o)^ p(m,o) used from p'"''^'. For example, as 

i:ri 

for the block shown in Fig. 9, comprising four pixels with 
py their pixel intensity values indicated inside, images p'"'~^'°', 

Q p(m-i,i)^ p(m-i,2) p{m-i,3) acquire "3", -8", -6" and "10", 

ru 

respectively, according to the rules described in [1.2]. This 
|>g 15 block at the m-th level is replaced at the (m-l)th level by 
respective single pixels thus acquired. Therefore, the size 
of the subimages at the (m-l)th level is 2"'"^X2'""\ 

After m is decremented (S103 in Fig. 8), it is ensured 
that m is not negative (S104). Thereafter, the process 
20 returns to SlOl, so that subimages of the next level of 

resolution, i.e., a next coarser level, are generated. The 
above process is repeated until subimages at m=0 (0-th level) 
are generated to complete the process at SIO. The size of the 
subimages at the 0-th level is IX 1. 
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Fig. 10 shows source hierarchical images generated at 
SIO in the case of n=3. The initial source image is the only 
image common to the four series followed. The four types of 
subimages are generated independently, depending on the type 
5 of critical point. Note that the process in Fig. 8 is common 
to Sll shown in Fig. 1, and that destination hierarchical 
images are generated through a similar procedure. Then, the 
process at SI in Fig. 6 is completed. 
J;* In this base technology, in order to proceed to S2 shown 

y 10 in Fig. 6 a matching evaluation is prepared. Fig. 11 shows 
j'p the preparation procedure. Referring to Fig. 11, a plurality 

py of evaluation equations are set (S30) . The evaluation 

[j^ equations may include the energy C^"''^ concerning a pixel 

i^p value, introduced in [1.3.2.1], and the energy D^f'^^ concerning 

i'U 15 the smoothness of the mapping introduced in [1.3.2.2]. Next, 
by combining these evaluation equations, a combined evaluation 
equation is set (S31) . Such a combined evaluation equation 
may be + Df-'K Using ti introduced in [1.3.2.2], we have 

zZ(^s?+^o;;:^^+^!s?) — (52) 

20 In the equation (52) the sum is taken for each i and j where i 
and j run through 0, 1,... , 2"" "\ Now, the preparation for 
matching evaluation is completed. 

Fig. 12 is a flowchart showing the details of the 
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process of S2 shown in Fig. 6. As described in [1], the 
source hierarchical images and destination hierarchical images 
are matched between images having the same level of 
resolution. In order to detect global correspondenge 
correctly, a matching is calculated in sequence from a coarse 
level to a fine level of resolution. Since the source and 
destination hierarchical images are generated using the 
critical point filter, the location and intensity of critical 
points are stored clearly even at a coarse level. Thus, the 
result of the global matching is superior to conventional 
methods . 

Referring to Fig. 12, a coefficient parameter r\ and a 
level parameter m are set to 0 (S20) . Then, a matching is 
computed between the four subimages at the m-th level of the 
source hierarchical images and those of the destination 
hierarchical images at the m-th level, so that four types of 
submappings f*"''^' (s=0, 1, 2, 3) which satisfy the BC and 
minimize the energy are obtained (S21) . The BC is checked by 
using the inherited quadrilateral described in [1.3.3]. In 
that case, the submappings at the m-th level are constrained 
by those at the (m-l)th level, as indicated by the equations 
(17) and (18) . Thus, the matching computed at a coarser level 
of resolution is used in subsequent calculation of a matching. 
This is called a vertical reference between different levels. 
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If m=0, there is no coarser level and this exceptional case 
will be described using Fig. 13. 

A horizontal reference within the same level is also 
performed. As indicated by the equation (20) in [1.3.3], 
f(m,3)^ f(m,2) ^(m,i) respectively determined so as to be 

analogous to ft"''^), fC^^D and ft-^-o'. This is because a 
situation in which the submappings are totally different seems 
unnatural even though the type of critical points differs so 
long as the critical points are originally included in the 
same source and destination images. As can been seen from the 
equation (20) , the closer the submappings are to each other, 
the smaller the energy becomes, so that the matching is then 
considered more satisfactory. 

As for f<™'°>^ which is to be initially determined, a 
coarser level by one may be referred to since there is no 
other submapping at the same level to be referred to as shown 
in the equation (19) . In this base technology, however, a 
procedure is adopted such that after the submappings were 
obtained up to f*"''^), f(m,o) recalculated once utilizing the 
thus obtained subamppings as a constraint. This procedure is 
equivalent to a process in which s=4 is substituted into the 
equation (20) and f'™'^* is set to f*™'"* anew. The above 
process is employed to avoid the tendency in which the degree 
of association between f^™'°> and f'"''^^ becomes too low. This 
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scheme actually produced a preferable result. In addition to 
this scheme, the submappings are shuffled in the experiment as 
described in [1.7.1], so as to closely maintain the degrees of 
association among submappings which are originally determined 
independently for each type of critical point. Furthermore, 
in order to prevent the tendency of being dependent on the 
starting point in the process, the location thereof is changed 
according to the value of s as described in [1.7]. 

Fig. 13 illustrates how the submapping is determined at 
the 0-th level. Since at the 0-th level each sub-image is 
consitituted by a single pixel, the four submappings f are 
automatically chosen as the identity mapping. Fig. 14 shows 
how the submappings are determined at the first level. At the 
first level, each of the sub-images is constituted of four 
pixels, which are indicated by solid lines. When a 
corresponding point (pixel) of the point (pixel) x in p'^'^' is 
searched within q'^'^', the following procedure is adopted: 

1. An upper left point a, an upper right point b, a lower left 
point c and a lower right point d with respect to the point x 
are obtained at the first level of resolution. 

2. Pixels to which the points a to d belong at a coarser level 
by one, i.e., the 0-th level, are searched. In Fig. 14, the 
points a to d belong to the pixels A to D, respectively. 
However, the pixels A to C are virtual pixels which do not 

MN-70022 



54 

exist in reality. 

3. The corresponding points A' to D' of the pixels A to D, 
which have already been defined at the 0-th level, are plotted 
in q*^'^*. The pixels A' to C are virtual pixels and regarded 
to be located at the same positions as the pixels A to C. 

4. The corresponding point a' to the point a in the pixel A is 
regarded as being located inside the pixel A' , and the point 
a' is plotted. Then, it is assumed that the position occupied 
by the point a in the pixel A (in this case, positioned at the 
lower right) is the same as the position occupied by the point 
a' in the pixel A' . 

5. The corresponding points h' to d' are plotted by using the 
same method as the above 4 so as to produce an inherited 
quadrilateral defined by the points a' to d' . 

6. The corresponding point x' of the point x is searched such 
that the energy becomes minimum in the inherited 
quadrilateral. Candidate corresponding points x' may be 
limited to the pixels, for instance, whose centers are 
included in the inherited quadrilateral. In the case shown in 
Fig. 14, the four pixels all become candidates. 

The above described is a procedure for determining the 
corresponding point of a given point x. The same processing 
is performed on all other points so as to determine the 
submappings. As the inherited quadrilateral is expected to 
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become deformed at the upper levels (higher than the second 
level) , the pixels A' to D' will be positioned apart from one 
another as shown in Fig. 3. 

Once the four submappings at the m-th level are 
determined in this manner, m is incremented (S22 in Fig. 12) . 
Then, when it is confirmed that m does not exceed n (S23) , 
return to S21. Thereafter, every time the process returns to 
S21, submappings at a finer level of resolution are obtained 
until the process finally returns to S21 at which time the 
mapping f^''' at the n-th level is determined. This mapping is 
denoted as f (^=0) because it has been determined relative to 
Ti=0. 

Next, to obtain the mapping with respect to other 
different r|, t| is shifted by Lrj and m is reset to zero (S24). 
After confirming that new ri does not exceed a predetermined 
search-stop value Timax{S25), the process returns to S21 and the 
mapping f^'" (7 = ^7) relative to the new ti is obtained. This 
process is repeated while obtaining f ''^^ (// = /A77) (/ = 0,1,...) at S21. 
When T| exceeds rimax, the process proceeds to S2 6 and the optimal 
■n='nopt is determined using a method described later, so as to 
let f ^"M'n=Tiopt) be the final mapping f'^'L 

Fig. 15 is a flowchart showing the details of the 
process of S21 shown in Fig. 12. According to this flowchart, 
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the submappings at the ra-th level are determined for a certain 
predetermined x\. In this base technology, when determining 
the mappings, the optimal X is defined independently for each 
submapping . 

Referring to Fig. 15, s and A, are first reset t.o zero 

(5210) . Then, obtained is the submapping f that minimizes 
the energy with respect to the then X (and, implicitly, t|) 

(5211) , and the thus obtained submapping is denoted as 

f (in,s) (^=0) . In order to obtain the mapping with respect to 
other different X, X is shifted by A/l . After confirming that 
the new A, does not exceed a predetermined search-stop value A^ax 

(5213) , the process returns to S211 and the mapping f ^'^'^^ 
(A, = AZ) relative to the new A, is obtained. This process is 
repeated while obtaining f^'^'^^(A. = iAA)(i = 0,l,-)- When X exceeds 
T^xr the process proceeds to S214 and the optimal X=Xopt is 
determined , so as to let f {X=Xopt) be the final mapping f f'^'^' 

(5214) . 

Next, in order to obtain other submappings at the same 
level, Xis reset to zero and s is incremented (S215) . After 

confirming that s does not exceed 4 (S216) , return to S211. 
When s=4, f*'^'^* is renewed utilizing f'"''^* as described above 
and a submapping at that level is determined. 

Fig. 16 shows the behavior of the energy C^"'-'^ 
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corresponding to f '""'^^ (A = zAA) (/ = 0, !,...)« for a certain m and s 
while varying X,. As described in [1.4], as ^.increases, Cf'^'^ 
normally decreases but changes to increase after X exceeds the 
optimal value. In this base technology, X in which Cj?""'^ 
5 becomes the minima is defined as Xopt- As observed in Fig. 16, 
even if Cf''^ begins to decrease again in the range A>Xopt , the 
mapping will not be as good. For this reason, it suffices to 
l;^' pay attention to the first occurring minima value. In this 

7^ base technology, X,opt is independently determined for each 

m 

„P 10 submapping including f'^'. 

f 

^■^ Fig. 17 shows the behavior of the energy Cf^ 

corresponding to f (;/ = r'A^) (/ = 0,1,...) while varying T]. Here 
too, Cf^ normally decreases as rj increases, but Cf^ changes to 
increase after Tj exceeds the optimal value. Thus, t| in which 
15 Cf^ becomes the minima is defined as rjopt- Fig, 17 can be 
considered as an enlarged graph around zero along the 
horizontal axis shown in Fig. 4. Once T|opt is determined, f^"^' 
can be finally determined. 

As described above, this base technology provides 
20 various merits. First, since there is no need to detect 

edges, problems in connection with the conventional techniques 
of the edge detection type are solved. Furthermore, prior 
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knowledge about objects included in an image is not 
necessitated, thus automatic detection of corresponding points 
is achieved. Using the critical point filter, it is possible 
to preserve intensity and locations of critical points even at 
5 a coarse level of resolution, thus being extremely 
advantageous when applied to object recognition, 
characteristic extraction, and image matching. As a result, 
it is possible to construct an image processing system which 
I''* significantly reduces manual labor. 

^i.^. 10 Some further extensions to or modifications of the 

i;n 

above-described base technology may be made as follows: 
ry ■ (1) Parameters are automatically determined when the matching 
Q is computed between the source and destination hierarchical 

m 

images in the base technology. This method can be applied not 
i|j 15 only to the calculation of the matching between the 

hierarchical images but also to computing the matching between 
two images in general. 

For instance, an energy Eo relative to a difference in 
the intensity of pixels and an energy Ei relative to a 
20 positional displacement of pixels between two images may be 
used as evaluation equations, and a linear sum of these 
equations, i.e., Etot=aEo+Ei, may be used as a combined 
evaluation equation. While paying attention to the 
neighborhood of the extrema in this combined evaluation 
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equation, a is automatically determined. Namely, mappings 
which minirnize Etot are obtained for various a's. Among such 
mappings, a at which Etot takes the minimum value is defined as 
an optimal parameter. The mapping corresponding to this 
parameter is finally regarded as the optimal mapping between 
the two images. 

Many other methods are available in the course of 
setting up evaluation equations. For instance, a term which 
becomes larger as the evaluation result becomes more 
favorable, such as 1/Ei and I/E2, may be employed. A combined 
evaluation equation is not necessarily a linear sum, but an n- 
powered sum (n=2, 1/2, -1, -2, etc.), a polynomial or an 
arbitrary function may be employed when appropriate. 

The system may employ a single parameter such as the 
above a, two parameters such as Tj and X as in the base 
technology, or more than two parameters. When there are more 
than three parameters used, they may be determined while 
changing one at a time. 

(2) In the base technology, a parameter is determined in a 
two-step process. That is, in such a manner that a point at 
which C)"'*^ takes the minima is detected after a mapping such 
that the value of the combined evaluation equation becomes 
minimum is determined. However, instead of this two-step 
processing, a parameter may be effectively determined, as the 
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case may be, in a manner such that the minimum value of a 
combined evaluation equation becomes minimum. In this case, 
aEo+pEi, for example, may be used as the combined evaluation 
equation, where a+p=l may be imposed as a constraint so as to 
equally treat each evaluation equation. The automatic 
determination of a parameter is effective when determining the 
parameter such that the energy becomes minimum. 

(3) In the base technology, four types of submappings 
related to four types of critical points are generated at each 
level of resolution. However, one, two, or three types among 
the four types may be selectively used. For instance, if 
there exists only one bright point in an image, generation of 
hierarchical images based solely on f'"''^) related to a maxima 
point can be effective to a certain degree. In this case, no 
other submapping is necessary at the same level, thus the 
amount of computation relative on s is effectively reduced. 

(4) In the base technology, as the level of resolution of an 
image advances by one through a critical point filter, the 
number of pixels becomes 1/4. However, it is possible to 
suppose that one block consists of 3X3 pixels and critical 
points are searched in this 3X3 block, then the number of 
pixels will be 1/9 as the level advances by one. 

(5) In the base technology, if the source and the 
destination images are color images, they would generally 
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first be converted to monochrome images, and the mappings then 
computed. The source color images may then be transformed by- 
using the mappings thus obtained. However, as an alternate 
method, the subraappings may be computed regarding each RGB 
component . 

Preferred Embodlmen-bs for Multivariate Space Processing 

A technique for multivariate space processing utilizing 
the above-described base technology will now be described. In 
the following, the term "variate" or "variable" is used to 
describe each of a plurality of variable quantities in a 
multivariate space. It is to be noted that a variate may also 
sometimes be called a parameter and these terms may generally 
be treated in the same manner as the term "dimension", as the 
case may be. Further, when a variate is referred to as a 
variate of an object, the variate may also be defined as an 
attribute of the object. 

Fig. 18 shows a structure of a multivariate space 
processing apparatus 10 according to an embodiment of the 
invention. This apparatus includes: a preprocessing unit 12 
which acquires multivariate data, which in a particular case 
is object data (OD) , and degenerates the object data into 
three variates; a conversion unit 14 which receives the 
degenerated data and then generates a first image and a second 
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image from the degenerated data; a matching processor 16 which 
performs a matching computation on the first image II and the 
second image 12; a corresponding point file storage unit 18 
which stores a corresponding point file F obtained as a result 
5 of the matching computation; an intermediate image generator 
20 which generates an intermediate image based on the first 
image II, second image 12 and corresponding point file F; and 
a display control unit 22 which performs any processing 
M necessary for displaying the intermediate image thus 

;J 10 generated. The conversion unit 14 also generates an authentic 

■ intentiediate image AIF at a predetermined point between the 
1=11 first image II and the second image 12, and outputs this 

q authentic intermediate image AIF to a comparator 24. The 

ru 

M intermediate image generator 20 also outputs a virtual 

Q 15 intermediate image VIF generated from the first image II and 
the second image 12 to the comparator 24. The comparator 24 
compares these two intermediate images, and stores a 
comparison result in a comparison result file storage unit 26. 
Preferably, the preprocessing unit 12 also stores the 
20 three variates selected by a user in a variate-selected memory 
unit 28 for use in subsequent processing. Moreover, the 
preprocessing unit 12 also preferably has access to comparison 
results for previous data by referring to the comparison 
result file storage unit 26. In this way, the preprocessing 
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unit 12 can display to a user previously selected variates and 
an indication of which previously selected variates provided a 
good comparison result. The user can utilize this information 
to aid in selection of the three variates to be used for the 
current process. 

Fig. 19 is a flowchart showing an overview of the 
processing in the multivariate space processing apparatus 10. 
The object data OD are first input to the preprocessing unit 
12 where the input object data OD are degenerated into three 
variates (S300) . Preferably, the preprocessing unit 12 
displays previously selected variates and some indication of 
the order of approximation between the authentic intermediate 
images AIF and the virtual intermediate images VIE for the 
previously selected variates, as described above. This may 
assist the user in selecting a certain variate or variates to 
be targeted in the visualization. This information is 
particularly effective in a case where there are many variates 
of the same or similar level of importance. 

Moreover, the preprocessing unit 12 preferably displays 
to the user previously selected variates for each type of 
object data OD. Thereby, the variates that the user was 
interested in in the past or the variates that the user chose 
for a predetermined purpose can be provided to the user to 
allow easier and more efficient selection of variates by the 
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user. 

Degeneration by the preprocessing unit 12 can be 
realized by various methods. As a simple example, the three 
variates the user selected, as they are, are selected and the 
remaining variates are ignored. As another method, in the 
course of reducing the variates, a conversion process can be 
performed on one or more remaining variates at an arbitrary 
timing. In a case where such a conversion is performed, the 
converted variates may also be stored in the variate-selected 
memory unit 28. 

Similar to a case where the degeneration method involves 
a certain dimension being simply cut off, an alternative may 
be to assign a specific value for that dimension. For 
example, in a case where a three-dimensional space is 
degenerated into a two-dimensional space, the height may be 
set to some constant. Here, it suffices that the degeneration 
is such that n dimensions are degenerated (reduced) to a 
dimension number which is less than n. For example, a three- 
dimensional space may be converted to an arbitrary plane 
ax-i-by+cz+d=0 or the three-dimensional space may be converted 
to a curved surface such as ax^+by+c=0 . In any event, 
reducing the number of variates is called degeneration. 

After the object data OD have been degenerated into 
three variates, the degenerated data are input to the 
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conversion unit 14. The conversion unit 14 then waits for the 
user to input which variate is to be regarded as a reference 
among the three variates. The determination of which variate 
is to be used as a reference may also be set as a default or 
the like. Once a variate which will serve as the reference 

(simply referred to as a reference variate hereinafter) is 
specified, the conversion unit 14 generates the first image II 
and the second image 12 by setting two predetermined values to 
this reference variate (S302) . Again, the predetermined values 
may be set by the user, may be default values, or may be 
otherwise determined. 

As an example, the reference variate may be time t and 
the remaining two variates may be x and y spatial co- 
ordinates. By examining these variates, a first image at time 
t=tl and a second image at time t=t2 can be defined. If an 
image matching is computed between these two images, an 
approximate behavior of the two variates x and y in the 
interval t=[tl, t2] can be determined. Thus, it suffices to 
compute a matching between the two-dimensional images and the 
computational load can be reduced by selecting a suitable 
algorithm. The "base technology" provides a suitable 
algorithm and contributes to producing a suitable 
visualization . 

The first image II and the second image 12 thus 
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generated are output to the matching processor 16. The 
matching processor 16 performs a matching computation where, 
for example, the computation may focus on critical points, as 
described in the base technology (S304) . The corresponding 
point file F generated as a result of the matching computation 
is stored in the corresponding point file storage unit 18. 

The intermediate image generator 20 receives the first 
image II, the second image 12, and the corresponding point 
file F and generates a virtual intermediate image VIF by 
performing an interpolation computation based thereon. An 
example of a possible interpolation computation is described 
in more detail in the base technology. The intermediate image 
thus generated is output to the display control unit 22. The 
display control until performs any necessary processing of the 
intermediate image and the processed intermediate image is 
output to a display device (5308) . The virtual intermediate 
image VIF thus generated is also .output to the comparator 24 
S308) . 

Fig. 20 shows more detail of S300 shown in Fig. 19. The 
preprocessing unit 12 first refers to the comparison result 
file storage unit 2 6 and the variate-selected memory unit 28 
to determine whether or not reference information exits 
(S300A) . Here, "reference information" is used as a generic 
term for previous comparison results and previously selected 
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variates. If any reference information exists (Y in S300A) , 
it is displayed (S300B) whereas, if reference information does 
not exist (N in S300A) , the display of the reference 
information (S300B)is skipped and the process waits for the 
user to specify selected variates (S300C) . 

After the user selects three variates by either 
confirming the displayed reference information ' or by 
specifying new variates (Y in S300C) , the preprocessing unit 
12 degenerates the object data OD into the selected three 
variates . 

We now consider a stock price as an example of how a 
complicated multivariate object can be converted into a simple 
model by selecting three variates. As will be understood, 
there are generally a very large number of parameters that go 
into determining a stock price. Initially, three variates are 
tentatively selected based on a rule of thumb or an arbitrary 
selection after a search algorithm which has enumerated a 
number of possible variates. Next, a first image and a second 
image are generated and a virtual intermediate image thereof 
is generated by matching and interpolation. Thereafter, this 
virtual intermediate image and the actual stock price (which 
represents the authentic intermediate image) are compared to 
determine the suitability of the selection of the three 
variates. Thus, if the virtual intermediate image is close to 
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the authentic intermediate image, it is concluded that the 
three variates have high importance. This can be used or 
further tested by determining if it is possible to analyze or 
predict the stock price at another time by using these three 
5 variates. If the virtual intermediate image and the authentic 
intermediate image are not close or if further testing is 
needed, the virtual intermediate image and the authentic 
intermediate image may continue to be compared while the three 
variates are being changed, the selection of the reference 

10 variate is being changed and/or the value (s) set for the 
reference variate is/are being changed. 

Fig. 21 shows more detail of S302 shown in Fig. 19. The 
conversion unit 14, which has received the three variates, 
first inquires of the user as to which variate is to be 

15 regarded as the reference variate (denoted as t 

hereinafter) (S302A) . After the user specifies the reference 
variate at S302A, the conversion unit 14 requests two specific 
values of the reference variate t from the user (S302B) . 
These two values are necessary for generating the first image 

20 II and the second image 12. As noted above, the two values 

may alternatively be predetermined or calculated values. After 
the conversion unit 14 acquires the user specified two values 
t=tl and t=t2 at S302B, the conversion unit 14 determines the 
first image II and second image 12 (S302C) . Since values of 



the remaining two variates may be determined by substituting 
specific values into the reference variate t, the first image 
II and the second image 12 represented by these two variates 
can be determined efficiently. 

The processes performed by the matching processor 16 at 
S304 of Fig. 19 and the intermediate image generator 20 at 
S306 of Fig. 19 are generally as described in the base ' 
technology above. 

Fig. 22 shows more detail of S308 shown in Fig. 19. 
After the intermediate image has been generated by the 
intermediate image generator 20 {S306 of Fig. 19), the display 
control unit 22 requires this intermediate image and may also 
receive the first image II and second image 12 depending on 
the application. The display control unit 22 then converts 
the intermediate image into a data format which can be 
displayed on the display unit (S308A) . Thereby, a two- 
dimensional image is realized for visualization. 

The intermediate image is also sfent to the comparator 24 
by the intermediate image generator 20 as a virtual 
intermediate image VIF and the comparator 24 compares the 
authentic intermediate image AIF and the virtual intermediate 
image VIF (S308B) . For example, the comparator 24 may compare 
the authentic intermediate image AIF and virtual intermediate 
image VIF at a midpoint t=(tl+t2)/2 of t=tl and t=t2. In this 
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case the conversion unit 14 may be arranged to substitute 
t=(tl+t2)/2 as the reference variate and determine the 
remaining two variates at that value so that a two-dimensional 
image can be generated. This image is generated based on the 
actual data and is, in this sense, an authentic intermediate 
image. Similarly, the intermediate image generator 20 may use 
the same reference variate value t=(tl+t2)/2 so that a virtual 
intermediate image at that reference variate value is 
generated through an interpolation computation. 

The comparator 24 takes a difference of the authentic 
intermediate image AIF and the virtual intermediate image and 
determines whether the result is desirable or not according to 
an amount of the difference. The comparison result is stored 
in the comparison result file storage unit 26 (S308C) . 
Alternatively, the comparator 24 may simply compute the 
difference and may store this computed difference in the 
comparison result file storage unit 26. In this case, the 
preprocessing unit 12 may determine whether or not the 
comparison result is desirable. 

Figs. 23, 24 and 25 show experimental visualization 
results according to the present embodiment. Here, data on 
atmospheric pressure in the vicinity of the center of a 
tornado are considered as the object data. Fig. 23 and Fig. 
24 show distributions of atmospheric pressure viewed from a 
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transverse direction at predetermined times t=tl and t=t2, 
respectively. Thus, Fig. 23 may be considered the first image 
II and Fig. 24 may be considered the second image 12 described 
above. Fig. 25 shows a virtual intermediate image at 
5 t=(tl+t2)/2 between the first image II and the second image 
12. Now, suppose that the z axis is placed in the 
longitudinal direction of the tornado, then the first image II 
and the second image 12, shown in Fig. 23 and Fig. 24, 
i;* respectively, can be thought of as three-variate made up of 

10 (1) two dimensional information acquired when the object data 

m 

„P are projected on a plane parallel to the z axis and (2) the 

py time. Thus, the processing of the preprocessing unit 12 in 

Q this example is one in which the object data are projected on 

M a plane parallel to the z axis and the time parameters are 

15 .stored. 

ru 

Further, the processing of the conversion unit 14 can be 
considered as selecting time t as the reference variate and 
obtaining two dimensional images at the times tl and t2 . As 
can be seen from the result shown in Fig. 25, visualization 
20 results can be obtained efficiently by utilizing a simplified 
method according to the present embodiment. 

In other words, the method and apparatus for 
multivariate space processing according to the above 
embodiment are implemented in the tornado example such that 
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(1) the first image and the second image are acquired by 
projecting three-dimensional data of an object on a 
predetermined x-y plane and (2) a matching between the first 
image and the second image thus acquired is computed so that 
intermediate states in the range of time t=[tl, t2] can be 
visualized by interpolating the two-dimensional images at time 
t=tl and time t=t2. Computationally, once the matching is 
completed, there will not be much time needed for the 
interpolation. Moreover, since the interpolation can be 
performed at any point by arbitrarily varying the time t, the 
computational load becomes extremely low compared to the 
conventional visualization where a three-dimensional 
computation is performed sequentially. Thus, the present 
embodiment is effective particularly when one wants to quickly 
visualize the result showing a simulation in as simplified and 
convenient a manner as possible. 

As is seen in the tornado example, there is an implicit 
relationship between degeneration and projection of object 
data. That is, when reducing or degenerating the parameters 
of a multidimensional phenomenon or object, the object may be 
^^viewed" from various "angles", an inspection may be carried 
out to determine an angle at which a most characteristic image 
is obtained, and the degeneration to view the object at that 
particular angle may be repeated. Using a simple example, 
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consider a case where a three-dimensional object is projected 
as a "shadow picture" on a plane. For most three-dimensional 
objects, there will be an angle at which the original object 
can be relatively easily grasped from the shadow picture. A 
projection operation using this angle corresponds to the 
degeneration . 

Though the present invention has been described based on 
the above embodiments, the present invention is not limited to 
these specific embodiments alone and various modifications 
thereto are also effective as encompassed by the present 
invention. Some of such the modifications will be described 
here . 

In the present embodiment, the tornado example was used 
to illustrate a phenomenon in three-dimensional space with the 
added dimension of time. However, the present invention is 
not limited to three or four dimensions and can also be 
applied to handle multivariate phenomena or objects. For 
example, in analyzing a complex economic phenomenon, a similar 
visualization process can be performed in a simple manner by 
initially specifying three variates that the user should pay 
attention to among various parameters. Then, by varying 
selection of the three variates, it is possible to determine 
which particular variates have a higher degre of importance in 
the economic phenomenon- Moreover, based on the comparison 

MN-70022 



74 

results of the comparator 24, the variate or variates which 
are important in the economic phenomenon can be more easily 
determined on a case-by-case basis. 

In other words, by implementing the multivariate space 
processing method according to the present invention, it is 
possible to avoid the complex modeling which is generally 
required in conventional simulation and analysis of 
complicated multivariate objects or phenomena, and 
visualization can be realized using a simple and convenient 
method. Thus, one of the main advantages of the embodiments 
of the invention lies in the proposition and realization of 
this new methodology by which arbitrary phenomena can be 
understood by a process of relatively simple matching of two- 
dimensional images and without modeling. 

Although the present invention has been described by way 
of exemplary embodiments, it should be understood that many 
changes and substitutions may be made by those skilled in the 
art without departing from the spirit and the scope of the 
present invention which is defined by the appended claims. 
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